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IMPLICIT APPROXIMATE-FACTORIZATION SCHEMES FOR THE EFFICIENT SOLUTION 


OF STEADY TRANSONIC FLOW PROBLEMS 

W. F, Ballhaus,* A. Jameson, ^ and J. Albert"^ 

Ames Research Center 
and 

Ames Directorate 

U.S. Army Air Mobility R&D Laboratory 


INTRODUCTION 


The objective of this work was to determine the feasibility of using 
implicit approximate-factorization algorithms (AF) , instead of successive 
line over-relaxation algorithms (SLOR) , to solve steady-state transonic flow 
problems. In this investigation, two AF solution procedures are compared 
with SLOR in terms of efficiency, reliability, and flexibility. 

Murman and Cole (ref. 1) developed the first computer-programmable SLOR 
algorithm for the solution of transonic flow problems. That solution proce- 
dure proved to be about an order of magnitude more efficient computationally 
than the first computer-programmed transonic flow algorithm, the explicit, 
time-accurate procedure of Magnus and Yoshihara (ref. 2). Since that time, 
the Murman-Cole procedure has been improved, extended, and applied to a 
variety of aerodynamic problems (reviews are given in refs. 3 and 4). SLOR 
algorithms have generally proved to be reliable, in the sense of convergence, 
and flexible, in the sense that they can be easily adapted to a wide range of 
applications . 

A practical limitation on the class of problems that can be treated is 
the computer time required to obtain a converged solution. Hence a number 
of potentially more efficient iterative solution procedures have been pro- 
posed as alternatives to SLOR. These are semidirect methods using fast 
Poisson solvers, extrapolation, the multigrid approach, and implicit approxi- 
mate factorizations. 

In recent years fast direct methods have been developed for solving 
matrix equations for the discrete Laplacian (refs. 5-7). For the transonic, 
small-disturbance potential equation 
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this suggests an iteration procedure of the form 
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such that a discrete Poisson equation is solved for each iteration n. 

Such a scheme has been proposed by Martin and Lomax (ref. 8) and by Martin 
(ref. 9). Schemes of this type produce substantial increases in computational 
efficiency, especially for cases with embedded supersonic regions that are 
small. However, their application has been limited because the most efficient 
Poisson solver algorithms require uniform grid spacing in all but one coordi- 
nate direction, and they require special, complex treatment of internal 
boundaries. A hybrid Poisson solver-SLOR scheme proposed by Jameson (ref. 10) 
avoids these difficulties for some special cases of two-dimensional flows by 
introducing a coordinate mapping in a simple domain. 


Extrapolation has proved useful for the acceleration of SLOR convergence 
in cases in which a dominant eigenvalue exists and can be Identified (refs. 11 
and 12). It has also been used, in a somewhat different form, to accelerate 
convergence in the fast Poisson solver approach (refs. 8 and 9), 

The multigrid method was first proposed by Federenko (ref. 13), developed 
by Brandt (ref. 14), and recently applied to the solution of transonic flow 
fields by South and Brandt (ref. 15). Preliminary results indicate a sub- 
stantial increase in efficiency over SLOR for computations on uniform meshes. 
However, the method is sensitive to the aspect ratio of the mesh cells, 
leading to complications in the treatment of nonuniform grids. 


Here we test approximate-factorization schemes to determine if they are 
capable of providing accelerated convergence and whether there are any 
restrictions limiting their application to transonic flow calculations. 


APPROXIMATE FACTORIZATIONS 


General Requirements 

We seek the solution to the difference equation 

L<J> = (AS + a )$ = 0 (3) 

T xx yy T 

where 5 XX and Syy are second central difference operators, 4 is the 
solution vector, and A is a constant. An iterative solution procedure to 
solve (3) can be written 
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NC n + cL<j> n = 0 


where C n is the correction vector <|> n - q>‘‘ computed for cycle n+1. We 

will call R = L4> n the tesidual at the nth iteration, and o is a parameter 
to be determined. 
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Let e"‘ = t|> n - <J> be the error after the nth 
and equations (3) and (4) can be combined to give 
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For the iteration procedure to converge, the modulus of every eigenvalue of 
M must be less than unity. 


The type of an iterative solution procedure is determined by the choice 
of the operator N. To ensure convergence, N should be chosen to resemble 
the operator L aB closely as possible and, furthermore, N -1 L<() n must be 
easily computable. If N is identical to L, for example, then M = I - aN L 
vanishes (for a - 1), and an exact solution to L<f> “ 0 is obtained in a 
single cycle. In the special case of the discrete Laplacian in a rectangular, 
uniform, Cartesian grid, this can be accomplished by the use of one of the 
fast Poisson algorithms. It is usually impractical to invert L, however, 
and one must use an iteration scheme in which N is an approximation to L. 
For example, for SLOR 
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where i,j are grid point Indices in the x,y directions, E^C-^j = C i _ 1 ., 
and a is the relaxation parameter. Note that N contains the ’ 

Syy operator in L but contains only the lower diagonal part of the 6^ 
operator. Hence, SLOR will require more than one cycle to converge. 

Actually, the required number of cycles will increase as the number of x 
grid points increases, because a grid point is only influenced by a single 
grid point to the right of it in the x direction during one cycle. 


The underlying idea of the approximate factorization methods is to 
construct N as a product of two or more factors 


N = N l N 2 . . . N k 


each of which is restricted to a form leading to equations that can be easily 
solved. For example, each factor may have a block tridiagonal structure. 

On the other hand, the added flexibility resulting from the use of several 
factors allows L to be more closely approximated by N. Moreover, 
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algorithms can be constructed along these lines in such a way that the final 
N operator is fully Implicit, so that every grid point is influenced by every 
other grid point during each cycle. 


AF Scheme 1 

Considering equation (3) in the case when A > 0 so that the equation is 
elliptic, an approximate factorization can be realized as 

(a - 6 ) (a - Ad )C n « accL<j> n (7) 

yy xx 

where a is a parameter to be chosen, and a is an over-relaxation factor. 
Multiplying out the two factors on the left, the operator N, which should 
now be an approximation to aL, has the form 


N «= aL - a 2 I - 
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The success of the method will depend on aL dominating the error terms 

-a 2 I - <S VV A S . The required equations can easily be solved by letting 
y y xx 

f n = (a - Ad )C n and solving directly for f the tridiagonal matrix equa- 
tion XX 

(a - 6 yy )f n = crocL4> n (8) 

followed by a direct solution for C n of the tridiagonal matrix equation 


(a - A6 )C n = f n (9) 

xx 

To estimate the rate of convergence for this scheme (which is a reformu- 
lation of the Peaceraan-Rachford scheme) note that equation (7) can be written 
in the form 

(a - 6 ) (a - Ad )e n+1 = oaL e 11 (10) 

yy xx 

Now, assuming periodic boundary conditions and a uniform Cartesian grid, let 


e n (x,y) = ^ G n (p,q)e ipx e iqy (11) 

p>q=i 

Because equation (10) is linear, we need consider only a single arbitrary 
Fourier component. Substituting in (10) and rearranging gives, in the case 
when a=2, 
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where n “ qAy and £ = pAx. For stability, that is, for the error to 
approach zero in the iteration procedure, 

I el < i (13) 

and this condition is satisfied for all E, and n eigenvalues in (12) . 


The number of iterations required to reduce the residual to a predeter- 
mined amount depends on |b|j that is, fewer iterations are required to 
achieve a specified degree of convergence as |b| decreases. The convergence 
rate, therefore, depends strongly on the choice of a. The choice 


Ay' 


(1 - cos n) 


(14) 


gives B B 0 for the particular eigenvalue n. Hence, for a problem with k 
grid points in the y-dlrection, corresponding to k eigenvalues, the solu- 
tion process will converge to zero error after k cycles. 


Precise estimation of the eigenvalues is generally not practical. 
Instead, a repeating sequence of a’s can be used with each element of the 
sequence chosen to maintain small values of |$| in a given range of 
eigenvalues (see ref. 16). Maximum and minimum values of a are estimated 
and then used to form a geometric sequence to cover the entire eigenvalue 
spectrum. The lowest eigenvalue is given approximately by n ~ Ay, which, 
from equation (14), gives a ^ - 1. For the high frequency error components, 
n ~ tt corresponding to ~ 4/Ay 2 . The sequence 


a 


k 



(2k-l)/2M 


k = 1 , 2 , . . . , M 


(15) 


is then used repetitively during the course of the computation. The sequence 
is adjusted in the following section to optimize convergence for particular 
transonic flow field computations. 

In solving equation (1) for transonic flow computations, scheme (7) is 
suitable for the treatment of the subsonic part of the flow, where 
A = (k - <(> x ) >0. To extend the scheme to treat flows with embedded super- 
sonic zones, the factorization is modified at supersonic points, where 
A = (k - <f> x ) <0, to the form 

(-A5 - S )aC n = accLd> n (16) 

xx yy T 
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where ^ xx Is the upwind second difference operator. Here, because of the 
use of upwind differencing, it is feasible to shift the term A 6 X3{ to the 
first factor while retaining an easily invertible form, and the second factor 
is essentially eliminated. On setting c » 1, the result is a transition to 
the fully implicit marching scheme of Murman and Cole (ref, 1) in the super- 
sonic zone. Transonic operators are used at the sonic and shock points to 
maintain conservation form (ref. 17). A term f3o x ( 6 X is a first-order, first 
upwind difference operator) is added to the first factor as a stabilizer with 
3 - Ax. At the shock point, D 6 X is added to the second factor for the first 
fow cycles to prevent instabilities associated with shock motion (e.g., see 
ref. 3 or 18). 


AF Scheme 2 

Implicit approximate-factorization schemes for the low-frequency tran- 
sonic equation 


K. <J> „ a (k - 4 1 ) ‘I’ + (17) 

l T xt r x T xx T yy 

where K 1 is a constant, were investigated by Ballhaus and Steger (ref. 18). 
For the conservative schemes tested, an Instability appeared whenever the 
following two conditions were satisfied; ( 1 ) the differencing was shifted 
from upwind to central across a shock and ( 2 ) the shock propagated at a rate 
greater than one spatial grid point per time step. However, a nonconservative- 
in-tlme (but conservative in the steady state) scheme was found to be much 
less susceptible to this instability. For this scheme it is not necessary 
to add D? x at the shock to maintain stability as for AF scheme 1. This is 
an attractive feature, since the appropriate value of D is difficult to 
determine. The scheme, a factorization of (17) which we will call AF 
scheme 2 (ref. 18), is 

[c. - (1 - - « yJ ,Ic n 

- {•*„ + <* x ia - VV* + ViVA 1 !*” a8) 

where ej “ [£] for Aj[£] zero, Aj = k - (<J>j +1 - < J ) j_ 1 )/ 2Ax « a a ( At ) _1 » 
and ~6 , X, are first-order, first forward, and backward difference operators. 

X X 

Note that the algorithm can be coded so that the <j) array need be stored for 
only a single time level (ref. 18). Note also that for subsonic flow 
(Aj >0), AF scheme 2 can be derived from AF scheme 1, equation (7), by 

replacing a by ot5 x and then removing <5 X from the second factor and the 
right-hand side. 

Substitution of form (11) Into (18) leads to an error growth factor 8 
that is not real. Choosing a to minimize 8 is therefore more difficult,, 
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and, in fact, there is no choice of a that can cause 3 to be zero for a 
given eigenvalue. Approximate values of ct£ and au for equation (15) are 
1 and Ay -1 . Of course, there is no guarantee that (15) is the optimal 
functional form for the acceleration parameter sequence. 

Thus, it appears at first glance that AF scheme 2 cannot be so highly 
tuned for rapid convergence as AF scheme 1, but it shows potential for greater 
convergence reliability because of its treatment of shock waves. 


OPTIMIZED ACCELERATION PARAMETERS 


Determination of the optimum relaxal on parameter for SLOR, or the 
optimum acceleration parameter sequence for approximate factorizations, for a 
practical computation is a formidable task to achieve analytically. However, 
comparisons of both optimum and non-optimum SLOR and approximate factorizations 
would be useful In assessing the relative merits of these solution procedures. 
Therefore, we have selected particular representative transonic flow problems 
for which "optimum” acceleration parameters can be determined by numerical 
optimization. 

The numerical optimization problem is formulated in the following way. 
First, we select as the objective , that is, the quantity to be minimized, 
a combination of parameters that represents the computational efficiency 
associated with a given set of decision variables: 

OBJ « (R n /Ri) l/N (19) 

OBJ is related to the average decrease in the residual per Iteration, R n is 
the maximum residual at the nth iteration, and N is the number of Iterations 
required to reduce the maximum residual to some specified value. 

Next, we define the decision variables, which are parameters that affect 
the convergence rate. For SLOR the decision variables are the subsonic 
relaxation parameters, one for each of the multiple meshes used to converge 
the solution. The supersonic relaxation parameter is always unity. For both 
AF schemes, the accleration parameter sequence is 

c* k = a^d/y) 2 ^ 1 k = 1,2,3, . . . , M (20) 

and the decision variable used in the optimization is y. The initial guess 

is y - (o^/^h) where and are given in the previous section. 

For coarse, medium, and fine grids, values of 4, 6, and 8, respectively, 
were used for M. 

The finite difference AF and SLOR codes were coupled with an executive 
optimization code, CONMIN (ref. 19). CONMIN is designed to minimize an 
objective, for a set of decision variables, subject to specified constraints. 
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Here the only constraint imposed was a limit of 2 on the relaxation parameter 
for SLOR. 

The present optimization formulation was designed to be as simple and 
economical as possible. Only single parameter optimizations are performed 
on each grid for both the AF and SLOR schemes. There is certainly no guaran- 
tee that either the functional form (20) or the arbitrarily chosen values of 
M used are optimum. 


COMPUTED RESULTS 


The convergence characteristics of the SLOR and AF schemes are investi- 
gated here for a simple test problem. The problem is that of computing the 
flow field about a nonlifting, parabolic-arc airfoil with a thickness-to- 
chord ratio of 10 percent. Three free-stream Mach numbers are considered: 

(1) Mo, “ 0.7, a subcrltlcal case; (2) M*, “ 0,84, a supercritical case; and 
(3) Moo = 0.9, a highly supercritical case. Computed surface pressure 
coefficients for these three cases are shown in figure 1. For all of the 
computations reported here, the solution process was considered completed 
when the absolute value of the maximum flow field residual dropped below 
10 -9 . Here the residual is defined as Ax 2 L<}>. 

The optimum convergence history of the AF and SLOR schemes for the 
Mo, = 0.7 case is shown in figure 2(a). This computation was performed 
using a single, fine grid with (128x32) (x,y) grid points. SLOR required 
350 iterations as opposed to 34 and 45 for AF schemes 1 and 2. The peaks in 
the AF schemes’ convergence histories correspond to the smaller values of 
ct(i.e., larger values of At) in the eight element sequence. 

SLOR convergence can usually be accelerated by the use of coarse-to-flne 
mesh interpolations. For example, converged solutions on a (64 X 16) grid 
were interpolated onto the fine (128x32) grid to provide a good starting 
solution. The resulting fine grid convergence history is illustrated in 
figure 2(b). Interpolation reduced the number of SLOR fine-mesh iterations 
from 350 to 141. Its effect on the AF schemes was small. Similar results 
are shown in figure 3 for the M ro = 0.84 case. The use of multiple meshes 
accelerated SLOR convergence here as in the M^ = 0.7 case. It also improved 
the AF-2 scheme convergence somewhat, but perhaps not enough to justify the 
additional complexity and computational work required in generating coarser 
grid solutions. However, the use of multiple meshes seems to be extremely 
important for AF-1. Good initial guesses approximately locate the shock wave 
and thus help prevent the "moving shock instability" mentioned in the 
preceding section. Difficulties with this type of instability were obtained 
in the "fine grid only" case, and hence results for that case <jr.e not 
reported in figure 3(a). For M m = 0.9, the supersonic bubble extends verti- 
cally from the airfoil surface to a point close to the far-field grid bound- 
ary. The optimum convergence history for this case is illustrated in fig- 
ure 4. Instabilities were encountered with AF-1 for both grid systems. In 
the fine grid only case, AF-2 required an order of magnitude fewer iterations 
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for convergence than SLOR. Coarse-to-f ine mesh interpolation cut the required 
number of SLOR iterations in half but had very little effect on the AF-2 
scheme convergence. The local peaks in the SLOR convergence history, fig- 
ure 4(b), correspond to changes in the location of the maximum residual. 

Up to this point we have compared the convergence performance of optimum 
SLOR and optimum AF, where here "optimum" refers to values of acceleration 
and relaxation parameters obtained from the single parameter numerical opti- 
mization procedure outline in the preceding section. However, in engineering 
work optimum parameters are almost never used. A good solution procedure 
must therefore provide fast, reliable convergence for a reasonable range of 
nonoptimum parameters. To test the nonoptimum behavior of the SLOR and AF 
schemes, we applied the optimum values obtained for the subsonic case, 

M m « 0.70, to the two supercritical cases, M^ « 0.84, 0.90; results are 
shown in figures 5 and 6 . Table 1 summarizes all of the results presented 
here in terms of number of Iterations required to achieve convergence. 

Table 2 gives optimum values of the acceleration and relaxation parameters. 
Results for coarse and medium grids are also Included. 

It is clear from Table 2 that the optimum acceleration parameter sequence 
is much more sensitive to the grid than to M„. For example, optimum values 
of y for AF-2 range from 1.363 to 1,455 for the fine grid over all three 
Mach numbers and both mesh systems. Hence, there is only slight deteriora- 
tion in convergence performance in the nonoptimum cases, and the effect of 
mesh system is very small. The nonoptimum SLOR results appear to be 
dominated by the difference between the optimum relaxation parameters for 
the fine grid only (ct 0 p T = 1.965) and for the fine grid with interpolation 
(oi 0 pT = 1.919). The Moo e 0.70 optimum acceleration parameter for the fine 
grid only case is nearly equal to that for Mm = 0.84. However, for the 
interpolated case, the optimum values for the two Mach numbers are substan- 
tially different. Hence, the fine grid only case converges faster, although 
less smoothly. Results for the AF-2 scheme at = 0.90 show the same trends 
as for M» = 0.84, that is, the convergence performance is (1) nearly equal to 
that for the optimum case and (2) insensitive to mesh interpolation. In the 
fine grid case, SLOR diverges because of the excessively large relaxation 
parameter . 


CONCLUDING REMARKS 


The present study indicates that implicit, approximate-factorization 
schemes can provide rapid and reliable convergence in finite-difference 
transonic flow computations. These schemes can be easily coded, require 
about the same storage as, and only about 50 percent more computational work 
per iteration than, present successive line over-relaxation algorithms. They 
can be used to obtain solutions on nonuniformly-spaced grids, a feature not 
presently available in some of the other convergence acceleration schemes. 

Of the two AF schemes investigated here, AF-1 converged the most 
rapidly. However, it was found to be very susceptible to the "moving shock 
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wave Instability" discussed in references 3 and 18. AF-2 was found to be 
insensitive to this instability for the cases computed here, and its conver- 
gence was found to be at least as reliable as that of SLOR and substantially 
faster. 

There should be no fundamental problems in applying the AF schemes to 
lifting cases. However, the application of these schemes to practical 
three-dimensional cases, or to the solution of the full potential equation, 
should be more difficult. The problem is to devise a stable, easily solvable, 
implicit factorization for the governing equations, which are complicated by 
the presence of cross derivatives. 
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TA3LE 1.- NUMBER OF ITERATIONS REQUIRED FOR CONVERGENCE 


M 

CO 

Grid 

system 

Parameters 

AF-1 

AF-2 

SLOR 

0.7 

Fine only 

Opt 

34 

45 

350 


Interpolated 

Opt 

29 

44 

141 

0.84 

Fine only 

Opt 

— 

84 

352 



Non Opt 

— 

99 

373 


Interpolated 

Opt 

51 

66 

253 



Non Opt 

58 

101 

508 

0.90 

Fine only 

Opt 

— 

82 

819 



Non Opt 

— 

91 

— 


Interpolated 

Opt 

— 

72 

411 



Non Opt 


93 

688 
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TABLE 2.- OPTIMUM ACCELERATION AND RELAXATION PARAMETERS 


M «- 

Crid 

system 

Y(AF-l) 

Y (AF-2) 

SLOR 

0.70 

Fine only 

1.813 


1.386 


1.965 


Interpolated 

(1) 2.276 

(1) 

1.576 

(1) 

1.830 



(2) 1.979 

(2) 

1.428 

(2) 

1.903 



(3) 1.843 

(3) 

1.363 

(3) 

1.919 

0.64 

Fine only 

— 


1.400 


1.961 


Interpolated 

(1) 2.319 

(1) 

1.730 

(1) 

1.821 



(2) 2.087 

(2) 

1.531 

(2) 

1.911 



(3) 1.895 

(3) 

1.455 

(3) 

1.960 

0.90 

Fine only 

— 


1.414 


1.926 


Interpolated 

(1) 

(1) 

1.765 

(1) 

1.858 



(2) 

(2) 

1.417 

(2) 

2 ,01 a 



(3) 

(3) 

1.404 

(3) 

1.976 


Noce: Numbers in parentheses correspond to mesh systems listed below: 


Mesh 

No. x points 

No. y points 

(1) 

32 

8 

(2) 

64 

16 

(3) 

128 

32 


Stability was maintained in some cases for relaxation parameters 
slightly greater than 2. 
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(b) Fine grid with mesh interpolation. 
Figure 2.- Concluded. 
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